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Abstract 

fi , Path integral Monte Carlo (PIMC) simulations have become an im- 

• 'T^ ' portant tool for the investigation of the statistical mechanics of quantum 

►^.,^1 systems. I discuss some of the history of applying the Monte Carlo method 

(-H , to non-relativistic quantum systems in path-integral representation. The 

O ,' principle feasibility of the method was well established by the early eight- 

ies, a number of algorithmic improvements have been introduced in the 

last two decades. 

T-H . 

>^ 

2 ■ 1 Introduction 

O. 

[~^ , Feynman's classic paper of 1948 presented a Space-Time Approach to Non- 

^^ ' Relativistic Quantum Mechanics, M or, in Hagen Kleinert's words, "an all- 

time global approach to the calculation of quantum mechanical amplitudes." 
Within the philosophy of this approach, we must find, as Kleinert often stressed, 
"^ , "all properties, including the Schrodinger wave functions, from the globally 

. ,^ ' determined time displacement amplitude." |p[ The Feynman path, governed by 

(— I ' study if we want to establish a truly independent alternative to Schrodinger's 

Oj. equation. In avoiding the operator formalism, the sum over paths provides an 

independent conceptual link between quantum and classical mechanics. 

Another attractive feature of the path integral formulation of quantum me- 
chanics is the bridge it allows to establish between quantum mechanics and sta- 
tistical mechanics. Technically, the oscillating exponential of the time-displace- 
ment amplitude turns into a positive Boltzmann weight if the paths are ex- 
pressed in imaginary time. The quantum mechanical propagator thus turns 
into a quantum statistical density matrix. It is this very feature which allows 

*To appear in: Fluctuating Paths and Fields. Festschrift Dedicated to Hagen Kleinert, 
Eds. W. Janke, A. Pelster, H.-J. Schmidt, and M. Bachmann (World Scientific, Singapore, 
2001). 
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the application of methods of classical statistical mechanics, notably the Monte 
Carlo method, to quantmii systems. 

The Monte Carlo method came into being roughly around the same time as 
the Feynman path. Anecdotally, the idea of gaining insight into a complex phe- 
nomenon by making various trials and studying the proportions of the respective 
outcomes occurred to Stanislaw Ulam while playing solitaire during an illness 
in 1946 P|. The immediate application was, of course, the problem of neutron 
diffusion studied in Los Alamos at that time. The name of the procedure first 
appeared in print in a classic paper by Metropolis and Ulam in 1949 Q, where 
the authors explicitly mentioned that the method they presented as a statis- 
tical approach to the study of integro-differential equations would sometimes 
be referred to as the Monte Carlo method. In classical statistical mechanics it 
quickly became a standard calculational tool. 

2 Where's Monte Carlo? 

The object of interest in Monte Carlo evaluations of Feynman's path integral is 
the quantum statistical partition function Z^ given, in operator language, as the 
trace of the density operator exp{—(3H) of the canonical ensemble (/3 = l/fceT') 
associated with a Hamilton operator describing N particles of mass rrii moving 
under the influence of a potential V, 

N ^2 
i—l 

Expressed as a Feynman integral, the density matrix elements read 

r{hl3)=r' / n/3 A 

(r|exp(-/3iJ)|r')= J I?r(r) exp J -i y L ({ri(T),rl(T)}) dr I (2) 

r(0)=r I ) 

where r = {fi, . . . , ttv}, L denotes the classical Lagrangian 

N 

L({ri(r),rl(r)}) = ^ I^^f ^ y(rl, . . . , f^(r)) (3) 

2=1 

expressed in imaginary time r.F] The particles are assumed to be distinguishable. 
To evaluate the trace, we only need to set r = r' and integrate over r. To take 
into account Bose or Fermi statistics for indistinguishable particles, the partition 
function splits into a sum of the direct Boltzmann part and parts with permuted 
endpoints. 

^ There have been attempts to apply the Monte Carlo method to path integrals also for 
real time. However, due to the oscillating exponential one then has to deal with problems of 
numerical cancellation, and it is much harder to obtain results of some numerical accuracy. 
Therefore, I shall here restrict myself to Monte Carlo work in imaginary time. 



The right hand side of Eq. (0) is a path integral for the 3N functions r. 
The idea of a Monte Carlo evaluation of this quantity is to sample these paths 
stochastically and to get (approximate) information about the quantum statis- 
tics of the system by averaging over the finite set of paths generated in the 
sampling process. 

Monte Carlo data always come with error bars and, in general, the errors 
associated with numerical Monte Carlo data stem from two distinct sources. A 
systematic error of Monte Carlo evaluations of the path integral follows from 
the need to identify the paths by a finite amount of computer information. 
This can be done by discretizing the paths at some set of points in the interval 
(0, ?i/3). For a single particle moving in one dimension, the simplest discrete 
time approximation for L time slices reads (e = UP/L) 



{x\ exp(-/?iJ)|a;') 
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where A = (2Trhe/m)^''^ and xq = x and xl = x' . Alternatively, one may 
expand the individual paths in terms of an orthogonal function basis, e.g. by 
the Fourier decomposition, 

, , (x' ~ x)t -^ . klTT 

and express the density matrix as 

(a;|exp(-/3ff)|a;') = lim Jexp-^ ^—{x - x')^) x 



L' 

n 



-^ / V{x{T))dT \ (6) 



where (Jk = [2,fi^ (3 / mijrk)'^]'^/'^ and J is the Jacobian of the transformation 
from the integral over all paths to the integral over all Fourier coefficients. A 
systematic error then arises from the loss of information by the finite number 
L of points Xi on the discretized time axis or by the finite number L' of Fourier 
components ak that are taken into account in the Monte Carlo sampling of the 
paths. 

The other error source of Monte Carlo data is the statistical error due to 
the finite number N^ of paths that form the sample used for evaluating the 
statistical averages. To make matters worse, the probability of configurations 
is, in general highly peaked, making an independent sampling of paths highly 
inefficient in most cases. The remedy is to introduce some way of "importance 
sampling" where configurations are generated according to their probability 
given by the exponential in Eq. (0). Statistical averages may then be computed 



as simple arithmetic means. A way to achieve this is by constructing Markov 
chains where transition probabiHties between configuration are constructed that 
aUow to generate a new configuration from a given one such that in the hmit 
of infinitely many configurations the correct probability distribution of paths 
results. A very simple and universally applicable algorithm to set up such a 
Markov chain is the Metropolis algorithm introduced in 1953 B. Here a new 
configuration is obtained by looking at some configuration with only one variable 
changed and accepting or rejecting it for the sample on the basis of a simple rule 
that depends only on the respective energies of the two configurations. The ad- 
vantages of importance sampling on the basis of Markov chains are obtained on 
the cost that, in general, successive configurations are not statistically indepen- 
dent but autocorrelated. The crucial quantity is the integrated autocorrelation 
time T^* of a quantity of interest O = (O) with O = (l/iV,„) Ej=i ^i and Oi 
computed for each path i in the sample. It enters the statistical error estimate 
Ao for expectation values of O computed from a Monte Carlo sample of Nm 
autocorrelated configurations as 




(7) 

where <7q , is the variance of Oi . 

With Monte Carlo generated samples of Feynman paths one can thus "mea- 
sure" thermodynamic properties of quantum systems like the internal energy 
and the specific heat, but also gain more detailed information about correlation 
functions, probability distributions and the like. In the low-temperature limit, 
P — > oo, quantum mechanical ground state properties are recovered. 

3 Blazing Trails 

A pioneer in the application of the Monte Carlo method to physics problems, 
notably by applying it to the Ising model, was Lloyd D. Fosdick. He ap- 
pears to have also been one of the first to consider the stochastic sampling 
of paths. In 1962, he considered the possibility of sampling paths |^ for what 
he called the conditional Wiener integral, i.e. the Wiener integral for fixed end 
points. As a toy example he investigated the expectation value of the func- 
tional exp — J„ L tt' x{T)x{T')dTdT' for a conditional Wiener process, and, 

more generally, for the quantity exp — J„ Vdr , i.e. he considered direct com- 
putation of the partition function for a quantum particle moving in a poten- 
tial V . He introduced a Fourier decomposition of the paths and generated 
these by direct Monte Carlo sampling of the Fourier components as Gaus- 
sian stochastic variables. He did some explicit sampling of his toy model to 
demonstrate the feasibility of the method but a theoretical consideration of 
the one-dimensional harmonic oscillator was not considered worthwhile to be 
put on the computer, even though Fosdick at the time was at the University 



of Illinois and had access to the university's ILLIAC computer. His examples 
were primarily used to investigate the principle feasibility and possible accu- 
racy obtainable by the method. Instead, Fosdick went along to consider a pair 
of two identical particles and presented some numerical results for this prob- 
lem. Continuation of the work on the two-particle problem together with a 
graduate student led to the publication of a paper on the Slater sum, i.e. the 
diagonal density matrix elements, for He^ in 1966 M, and on three-particle ef- 
fects in the pair distribution function for He'' in 1968 g. In the same year 
[^, Fosdick elaborated on the numerical method in a report on the Monte 
Carlo method in quantum statistics in the SIAM review. Instead of sampling 
the Fourier components he now used the discrete time approximation of the 
paths. Sampling of Xi at the discrete points was done using a trick that later 
would gain prominence in PIMC simulations in an algorithm called staging. The 
idea is to express the discretized kinetic term in the relative probablity density 
p{xi\xi-i,Xi+i) = (l/27re)exp [-{xt - Xi-i)^/2e - {xi+i - Xi)^/2e] as -{xi - 
Xi)'^/2a'^ + {xi-i —Xi+i)'^/Ae with Xi = (x^-i +Xi^i)/2 and to sample {xi—Xi)/a 
as an independent Gaussian variable. The procedure could be iterated recur- 
sively for all Xi and thus allowed to obtain statistically independent paths which 
were used to "measure" the potential energy term exp — Jq V{x(T))dT . 

In 1969, Lawande, Jensen, and Sahlin introduced Metropolis sampling of 
the paths in discrete time, broken line approximation [QOl. They investigated 
the ground state wave functions of simple one-dimensional problems (harmonic 
oscillator, square well, and Morse potential) and, theoretically, also addressed 
the problem of extracting information about excited energies and of simulating 
many particle problems. In a follow-up paper ||ll| they presented investigations 
of the Coulomb problem using Monte Carlo simulations of the path integral 
in polar coordinates. Not surprisingly, the singularity at the origin had to be 
avoided by artificial constraints and the authors admitted that a more rigorous 
justification of their procedure was called for. The path integral was later solved 
exactly by Duru and Kleinert in 1979 [^. It became clear that there were 
fundamental problems with such singularities in time-sliced path integrals g] . 

Little activity is recorded in the seventies, and I am only aware of a brief 
theoretical consideration of the possibility of Monte Carlo sampling of paths 
in a paper by Morita on the solution of the Bloch equation for many particle 
systems in terms of path integrals from 1973 ||l^. The paper is cited in a 
later one by J. A. Barker published in 1979 jl4| in which the one-dimensional 
particle in a box is considered, and numerical estimations of the ground state 
energy and wave function are presented. The data were obtained by introducing 
image sources to take account of the boundary conditions of the box and using 
Metropolis sampling of the broken line approximation of the paths. Incidentally, 
the analytic solution of this problem, i.e. of the path integral for the particle 
in the box was given by Janke and Kleinert almost simultaneously llq] . Barker 
also computed distribution functions for the problem of two particles in a box. 

Very much in the spirit of Lawande et al. but possibly unaware of this work, 
Creutz and Freedman published a didactic paper on a statistical approach to 



quantum mechanics in 1981 Q. They, too, performed Metropolis sampling of 
paths in the broken line approximation and studied the energies and ground 
state wave functions of the one-dimensional harmonic oscillator. The back- 
ground of these authors were Monte Carlo studies of gauge field theories, and 
the paper was meant as an attempt to better understand the Monte Carlo 
method by applying it to simple one-degree-of-freedom Schrodinger problems. 
It still is a useful introduction to the basics of the technique, and in particular it 
presents a brief primer on the theory of Markov chains underlying the Metropo- 
lis algorithm. To compute energies they introduced an alternative estimator by 
invoking the virial theorem. They also studied double well problems, presenting 
snap shot pictures of double-kink instanton paths. The problem of determining 
the energy level splitting was addressed by computing correlation functions. 

The papers by Lawande et al. and by Creutz and Freedman appear to have 
been cited very rarely, possibly because they presented their work as being only 
of pedagogic value and not so much because the Monte Carlo method could be 
a useful method to obtain numerical results for Schrodinger problems which, in 
real life, should be handled by numerical methods more suitable in this simple 
case. These remarks also hold for work published a little later by Shuryak 

Fosdick's work from 1962 was done very much at the forefront of the techno- 
logical possibilities of high-speed computing at the time. By the mid-eighties, 
path integral simulations of simple quantum mechanical problems had become 
both conceptually and technically "easy." Indeed, the exposition by Creutz 
and Freedman was already written in an introductory, didactic manner, and in 
1985 the simulation of the one-particle harmonic oscillator was explicitly pro- 
posed as an undergraduate project, to be handled on a Commodore CBM3032 
microcomputer, in a paper published in the American Journal of Physics p9| . 

4 Speeding up 

The feasibility of evaluating the quantum statistical partition function of many- 
particle systems by Monte Carlo sampling of paths was well established by the 
early eighties and the method began to be applied to concrete problems, in 
particular in the chemical physics literature. It had also become clear that the 
method had severe restrictions if numerical accuracy was called for. In addi- 
tion to the statistical error inherent to the Monte Carlo method, a systematic 
error was unavoidably introduced by the necessary discretization of the paths. 
Attempts to improve the accuracy by algorithmic improvements to reduce both 
the systematic and the statistical errors were reported in subsequent years. The 
literature is abundant and rather than trying to review the field I shall only 
indicate some pertinent paths of development. 

In Fourier PIMC methods, introduced in 1983 in the chemical physics context 
by Doll and Freeman EG, pl|, the systematic error arises from the fact that 
only a finite number of Fourier components are taken into account. Here the 
systematic error could be reduced by the method of partial averaging [E2l E3| . 



In discrete time approximations arising from the short-time propagator or, 
equivalently, the high-temperature Green's function various attempts have been 
made to find more rapidly converging formulations. Among these are attempts 
to include higher terms in an expansion of the Wigner-Kirkwood type, i.e. an 
expansion in terms of h /2m. Taking into account the first term of such an 
expansion would imply to replace the potential term eV{xj-i) in (^) by |^, p5| . 



eV{xj^i) 



dyV{y). 



(8) 



This improves the convergence of the density matrix (^ (from even less than 
0[1/L)) Hj] to 0{1/L'^). For the full partition function, the convergence of 
the simple discretization scheme is already of order ©(l/L^) since due to the 
cyclic property of the trace, the discretization eV{xj^i) is then equivalent to a 
symmetrized potential term e{y{xj-i) + V{xjy)/2. The convergence behaviour 
of these formulations follows from the Trotter decomposition formula. 
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valid for non-commuting operators A and i? in a Banach space 1 27 , identifying A 

with the kinetic energy jS^p^ /2mi and B with the potential energy (3V{{xi}). 
More rapidly converging discretization schemes were investigated on the basis of 
higher-order decompositions. Unfortunately, a direct, "fractal" decomposition 
fESl of the form 



<A+B) 



lim 



=/3lfp"2f„/32f 
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(10) 



inevitably leads to negative coefficients for higher decompositions [2^ and is 
therefore not amenable to Monte Carlo sampling of paths [^. Higher-order 
Trotter decomposition schemes involving commutators have proven to be more 
successful ipil p2l pi, M. In particular, a decomposition of the form 



lim Tr 



[IB.A\.B] 



(11) 



derivable by making use of the cyclic property of the trace, is convergent of 
order ©(l/L"*) and amounts to simply replacing the potential eV in (H) by an 
effective potential p2| 



''--''+IS<'")^- 



(12) 



Another problem for the numerical accuracy of PIMC simulations arises from 
the analog of the "critical slowing down" problem well-known for local update 
algorithms at second-order phase transitions in the simulation of spin systems 



and lattice field theory. Since the correlations (xjXj+k) between variables Xj 
and Xj+k in the discrete time approximation only depend on the temperature 
and on the gaps between the energy levels and not, or at least not appreciably, 
on the discretization parameter e, the correlation length C along the discretized 
time axis always diverges linearly with L when measured in units of the lattice 
spacing e. Hence in the continuum limit of e ^ with P fixed or, equivalently, 
of L — > CO for local, importance sampling update algorithms, like the standard 
Metropolis algorithm, a slowing down occurs because paths generated in the 
Monte Carlo process become highly correlated. Since for simulations using the 
Metropolis algorithm autocorrelation times diverge as |Q Tq^ oc L^ with z ~ 2 
the computational effort (CPU time) to achieve comparable numerical accuracy 
in the continuum limit L ^ oo diverges as L x L^ ^ i^+^. 

To overcome this drawback, ad hoc algorithmic modifications like introduc- 
ing collective moves of the path as a whole between local Metropolis updates 
were introduced then and again. One of the earliest more systematic and suc- 
cessful attempts to reduce autocorrelations between successive path configu- 
rations was introduced in 1984 by Pollock and Ceperly |Q. Rewriting the 
discretized path integral, their method essentially amounts to a recursive trans- 
formation of the variables Xi in such a way that the kinetic part of the energy can 
be taken care of by sampling direct Gaussian random variables and a Metropolis 
choice is made for the potential part. The recursive transformation can be done 
between some fixed points of the discretized paths, and the method has been 
applied in such a way that successively finer discretizations of the path were 
introduced between neighbouring points. Invoking the polymer analog of the 
discretized path this method was christened the "staging" algorithm by Sprik, 
Klein, and Chandler in 1985 fS^. 

The staging algorithm decorrelates successive paths very effectively because 
the whole staging section of the path is essentially sampled independently. In 
1993, another explicitly non-local update was applied to PIMC simulations 
[ p5[ p8| by transferring the so-called multigrid method known from the sim- 
ulation of spin systems. Originating in the theory of numerical solutions of 
partial differential equations, the idea of the multigrid method is to introduce a 
hierarchy of successively coarser grids in order to take into account long wave- 
length fluctuations more effectively. Moving variables of the coarser grids then 
amounts to a collective move of neighbouring variables of the finer grids, and 
the formulation allows to give a recursive description of how to cycle most effec- 
tively through the various levels of the multigrid. Particularly successful is the 
so-called W-cycle. Both the staging algorithm and the multigrid W-cycle have 
been shown to beat the slowing down problem in the continuum limit completely 
by reducing the exponent z to z w p9| . 

Another cause of severe correlations between paths arises if the probabil- 
ity density of configurations is sharply peaked with high maxima separated by 
regions of very low probability density. In the statistical mechanics of spin sys- 
tems this is the case at a first-order phase transition. In PIMC simulations 
the problem arises for tunneling situations like, e.g., for a double well potential 
with a high potential barrier between the two wells. In these cases, an unbiased 



probing of the configuration space becomes difficult because the system tends 
to get stuck around one of the probability maxima. A remedy to this problem 
is to simulate an auxiliary distribution that is flat between the maxima and 
to recover the correct Boltzmann distribution by an appropriate reweighting of 
the sample. The procedure is known under the name of umbrella sampling or 
multicanonical sampling. It was shown to reduce autocorrelations for PIMC 
simulations of a single particle in a one-dimensional double well, and it can also 
be combined with multigrid acceleration Q . 

The statistical error associated with a Monte Carlo estimate of an observ- 
able O cannot only be reduced by reducing autocorrelation times r^''. If the 
observable can be measured with two different estimators Ui that yield the same 
mean [/> ' — (Ui) with O — limL_,oo U^ , the estimator with the smaller vari- 
ance tJ^. is to be preferred. Straighforward differentiation of the discretized 
path integral (|4|) leads to an estimator of the energy that explicitly measures 
the kinetic and potential parts of the energy by 

The variance of this so-called "kinetic" energy estimator diverges with L. An- 
other estimator can be derived by invoking the path analog of the virial theorem 




= -{x.V'ix,)}, (14) 



and the variance of the "virial" estimator 

L L 

1=1 i=l 

does not depend on L. In the early eighties, investigations of the "kinetic" 
and the "virial" estimators focussed on their variances Ell E3, p2|. Some years 
later, it was pointed out Gsl that a correct assessment of the accuracy also 
has to take into account the autocorrelations, and it was demonstrated that for 
a standard Metropolis simulation of the harmonic oscillator the allegedly less 
successful "kinetic" estimator gave smaller errors than the "virial" estimator. 
In 1989 it was shown Q] that conclusions about the accuracy also depend on 
the particular Monte Carlo update algorithm at hand since modifications of the 
update scheme such as inclusion of collective moves of the whole path affect the 
autocorrelations of the two estimators in a different way. A careful comparison 
of the two estimators which disentangles the various factors involved was given 
only quite recently |Q. Here it was also shown that a further reduction of 
the error may be achieved by a proper combination of both estimators without 
extra cost. 



5 Concluding Remarks 

Application of the Monte Carlo method to quantum systems is not restricted to 
direct sampling of Feynman paths but this method has attractive features. It is 
not only conceptually suggestive but also allows for algorithmic improvements 
that help to make the method useful even when the problems at hand requires 
considerable numerical accuracy. However, algorithmic improvements like the 
ones alluded to above have tended to be proposed and tested mainly for simple, 
one-particle systems. On the other hand, the power of the Monte Carlo method 
is, of course, most welcome in those cases where analytical methods fail. For 
more complicated systems, however, evaluation of the algorithms and control 
of numerical accuracy is also more difficult. Only recently, a comparison of the 
efficiency of Fourier- and discrete time-path integral Monte Carlo for a cluster of 
22 hydrogen molecules was presented |Q — and debated |47| ^. Nevertheless, 
path integral Monte Carlo simulations have become an essential tool for the 
treatment of strongly interacting quantum systems, like, e.g., the theory of 
condensed hehum M. 
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